setwd("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS")
sel1=list.files(pattern = ".snp")

i=1
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file$tumor_var_freq=as.numeric(gsub("%","",file$tumor_var_freq))/100
file$normal_var_freq=as.numeric(gsub("%","",file$normal_var_freq))/100
file=data.frame(file$unitID,file$normal_var_freq,file$tumor_var_freq)
str=unlist(strsplit(gsub(".snp","",sel1[i]),"_"))
names(file)=c("unitID",paste0(str[1],".var.freq"),paste0(str[2],".var.freq"))
rt=file

for(i in 2:length(sel1)){
  file=read.table(sel1[i],head=T,sep="\t")
  file=file[file$exits=="T",]
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file$tumor_var_freq=as.numeric(gsub("%","",file$tumor_var_freq))/100
  file$normal_var_freq=as.numeric(gsub("%","",file$normal_var_freq))/100
  file=data.frame(file$unitID,file$normal_var_freq,file$tumor_var_freq)
  str=unlist(strsplit(gsub(".snp","",sel1[i]),"_"))
  names(file)=c("unitID",paste0(str[1],".var.freq"),paste0(str[2],".var.freq"))
  rt=merge(rt,file,by="unitID",all = T)
}
write.table(rt,"./ASM_log/869727.all.snp.vaf",quote=F,row.names = F,sep="\t")
rt$up.num=rowSums(rt[,2:35]>=0.5,na.rm = T)
rt$down.num=rowSums(rt[,2:35]<0.5,na.rm = T)
rt$mean.vaf=rowMeans(rt[,2:35],na.rm = T)
rt$group.up.down=ifelse(rt$mean.vaf>0.5,"up","down")
rt1=data.frame(rt$unitID,rt[,36:39])
write.table(rt1,"./ASM_log/869727.all.snp.vaf.up.down",quote=F,row.names = F,sep="\t")

################################################修改如下：20210321
setwd("E:/1.甲基化分析/ASM/ASM_snp-onlyWGS")
rt=read.table("./ASM_log/869727.all.snp.vaf",head=T,sep="\t")

affect=c("X1T","M7","M5","M1","M47","M51","M49","M43","M41","M28","M27","M30","M29","M26","M25","M35","M36")
sel1=paste(affect,".var.freq",sep="")

unaffect=c("X2B","M8","M6","M2","M48","M2","M50","M44","M42","M18","M17","M20","M19","M22","M21","M40","M39")
sel2=paste(unaffect,".var.freq",sep="")
rt$affect.vaf=rowMeans(rt[,sel1],na.rm = T)
rt$unaffect.vaf=rowMeans(rt[,sel2],na.rm = T)
rt$chazhi=rt$affect.vaf-rt$unaffect.vaf

rt$ASM.alt.group=ifelse(rt$chazhi>0,"up","down")
rt1=data.frame(unitID=rt$unitID,rt[,36:39])
write.csv(rt1,"./ASM_log/869727.all.snp.vaf.up.down.20210321.csv",quote=F,row.names = F)
